class: center, middle, inverse, title-slide # R-Vorkurs WS 22/23 ## Teil 3 ### Jens Klenke ### 07.10.2022 --- class: inverse # Übersicht <br></br> 1. Verteilungen und Zufallszahlen 2. Grafiken 3. Lineare Regression --- ## Verteilungen und Zufallszahlen ### Verteilungen R kennt standardmäßig viele diskrete Verteilungen, z.B. .font80[ - Binomialverteilung (`binom`) - Geometrische Verteilung (`geom`) - Hypergeometrische Verteilung (`hyper`) - Poisson-Verteilung (`pois`) - ... ] -- .font80[ ...sowie stetige Verteilungen: - Normalverteilung (`norm`) - t-Verteilung (`t`) - Chi-Quadrat-Verteilung (`chisq`) - F-Verteilung (`f`) - Gleichverteilung (`unif`) - ... ] --- ## Verteilungen ### <svg viewBox="0 0 581 512" style="height:1em;position:relative;display:inline-block;top:.1em;fill:#004c93;" xmlns="http://www.w3.org/2000/svg"> <path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"></path></svg>-Funktionen .font90[ Für jede Verteilung stehen vier Funktionen zur Verfügung: - Berechnung von Punktwahrscheinlichkeiten/-dichten `\(\Rightarrow\)` Präfix .blue[d] (.blue[d]ensity) - Werte der Verteilungsfunktion `\(\Rightarrow\)` Präfix .green[p] (.green[p]robability) - Berechnen von Quantilen `\(\Rightarrow\)` Präfix .red[q] (.red[q]uantile) - Erzeugen von 'Zufallszahlen' `\(\Rightarrow\)` Präfix .orange[r] (.orange[r]andom) ] -- `\(\Longrightarrow\)` Name der Funktion = Präfix + Verteilungskürzel (siehe vorherige Folie) -- .font90[ #### Beispiel: - Quantil der Normalverteilung: .red[q] + norm = `qnorm()` - Verteilungsfunktion der Binomialverteilung: .green[p] + binom = `pbinom()` - etc. ] --- ## Verteilungen ### Beispiel &mdash; Chi-Quadrat-Verteilung <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-1-1.png" width="500" style="display: block; margin: auto;" /> -- ```r dchisq(12, df = 10) # Dichte an der Stelle x=12 ``` ``` ## [1] 0.06692631 ``` --- ## Verteilungen ### Beispiel &mdash; Chi-Quadrat-Verteilung .font90[ Analog erhalten wir für die gleiche Verteilung Werte der Verteilungsfunktion: ```r pchisq(8, df = 10) # Wert der Verteilungsfunktion an der Stelle x=8 ``` ``` ## [1] 0.3711631 ``` ] -- .font90[ ...ein bestimmtes Quantil ```r qchisq(0.5, df = 10) # 0.5-Quantil der Verteilung ``` ``` ## [1] 9.341818 ``` ] -- .font90[ ...sowie Zufallszahlen ```r rchisq(5, df = 10) # 5 chi^2-verteilte Zufallszahlen ``` ``` ## [1] 5.472688 13.539234 4.101045 5.381317 8.123000 ``` ] --- ## Verteilungen ### Seeds .font80[ Wiederholen wir den Code von Slide 6, erhalten wir andere zufällige Zahlen. ```r rchisq(5, df = 10) ``` ``` ## [1] 3.303185 11.972862 9.340643 16.874128 4.991263 ``` ] -- .font80[ Manchmal wollen wir Ergebnisse reproduzierbar machen. Dann muss ein sog. *seed* gesetzt werden. Damit können an jedem Computer die selben Zufallszahlen erzeugt werden.] -- .font80[ ```r set.seed(385) # 385 ist ein Beispiel -- probiert andere seeds aus! rchisq(5, df = 10) ``` ``` ## [1] 13.158928 14.475007 7.448275 18.419976 8.822344 ``` ```r set.seed(385) rchisq(5, df = 10) ``` ``` ## [1] 13.158928 14.475007 7.448275 18.419976 8.822344 ``` ] --- ## Verteilungen ### <code>sample()</code> .font80[ Neben der Erzeugung von Zufallszahlen aus vordefinierten Verteilungen, kann man auch 'eigene' Zufallsexperimente durchführen ] -- .font80[ #### Beispiel Münzwurf: ```r # Mögliche Ergebnisse definieren K_Z <- c('Kopf', 'Zahl') ``` ] -- .font80[ Fünfmaliges Werfen einer *fairen* Münze: ```r # (*standardwert* prob = 1/length(K_Z) ) sample(K_Z, size = 5, replace = TRUE, prob = c(0.5, 0.5)) ``` ``` ## [1] "Kopf" "Zahl" "Zahl" "Zahl" "Kopf" ``` ] -- .font80[ Fünfmaliges Werfen einer *unfairen* Münze: ```r sample(K_Z, size = 5, replace = TRUE, prob = c(0.8, 0.2)) ``` ``` ## [1] "Kopf" "Kopf" "Kopf" "Zahl" "Kopf" ``` ] --- class: exercise_slide ## Verteilungen ### Übungsaufgaben .font90[ 1. Sei `\(X\sim t(5)\)`. Berechnen Sie `\(P(X<6)\)`, `\(P(3<X\leq7)\)` und `\(P(X>4)\)`. 2. Berechnen Sie das 0.95-Quantil einer `\(F(4,5)\)`-verteilten Zufallsvariable. 3. Berechnen Sie die Wahrscheinlichkeit dafür, den Jackpot im Lotto zu gewinnen (d.h. 6 Richtige aus 49). Vernachlässigen Sie bei Ihrer Berechnung Zusatz- oder Superzahlen. - *Hinweis*: Benutze die hypergeometrische Verteilung. 4. Erzeugen Sie 20 `\(\chi^2(5)\)`-verteilte Zufallszahlen ohne (!) dabei die `rchisq()`-Funktion zu benutzen. *Hinweis*: `\(\chi^2(n)=\sum_{i=1}^n Z_i^2\)` mit `\(Z_i\sim\mathcal{N}(0,1)\)` für alle `\(i=1,...,n\)`. 5. Ziehen Sie 10 mal `\(n = 10000\)` standardnormalverteilte Zufallszahlen und erechnen Sie für jeden Durchlauf das arithmetische Mittel. Schauen Sie sich danach alle 10 Mittelwerte an. Was fällt Ihnen auf? Sind Sie überrascht? - **Zusatzaufgabe:** Führen Sie dieselbe Simulationsstudie mit Cauchy-verteilten Zufallszahlen (`rcauchy()`) durch. Was fällt Ihnen nun auf? Können Sie sich das Ergebnis erklären? 6. Zeigen Sie, dass das Integral über die Dichtefunktion einer `\(\chi^2(15)\)`-verteilten Zufallsvariable 1 ist. - *Hinweis:* Nutze `integrate()` zum Integrieren einer <svg viewBox="0 0 581 512" style="height:1em;position:relative;display:inline-block;top:.1em;fill:black;" xmlns="http://www.w3.org/2000/svg"> <path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"></path></svg>-Funktion ] --- ## Grafiken ### Basic Plot Mit <svg viewBox="0 0 581 512" style="height:1em;position:relative;display:inline-block;top:.1em;fill:black;" xmlns="http://www.w3.org/2000/svg"> <path d="M581 226.6C581 119.1 450.9 32 290.5 32S0 119.1 0 226.6C0 322.4 103.3 402 239.4 418.1V480h99.1v-61.5c24.3-2.7 47.6-7.4 69.4-13.9L448 480h112l-67.4-113.7c54.5-35.4 88.4-84.9 88.4-139.7zm-466.8 14.5c0-73.5 98.9-133 220.8-133s211.9 40.7 211.9 133c0 50.1-26.5 85-70.3 106.4-2.4-1.6-4.7-2.9-6.4-3.7-10.2-5.2-27.8-10.5-27.8-10.5s86.6-6.4 86.6-92.7-90.6-87.9-90.6-87.9h-199V361c-74.1-21.5-125.2-67.1-125.2-119.9zm225.1 38.3v-55.6c57.8 0 87.8-6.8 87.8 27.3 0 36.5-38.2 28.3-87.8 28.3zm-.9 72.5H365c10.8 0 18.9 11.7 24 19.2-16.1 1.9-33 2.8-50.6 2.9v-22.1z"></path></svg> ist es einfach Grafiken zu erstellen, z. B. einen *Scatterplot*... -- ```r plot(rnorm(10)) ``` <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-12-1.png" width="600" style="display: block; margin: auto;" /> --- ## Grafiken ### Basic Plot ...oder einen Linienplot in roter Farbe -- ```r plot(rnorm(10), type = 'l', col = 'red') ``` <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-13-1.png" width="600" style="display: block; margin: auto;" /> --- ## Grafiken ### Cosmetics Neben Argumenten wie `col` können weitere Eigenschaften des Plots angepasst werden, z.B.: -- - `xlab`/`ylab`: Beschriftung von x-/y-Achse - `xlim`/`ylim`: Darzustellender Bereich von x-/y-Achse - `main`/`sub`: Titel/Untertitel der Grafik - `pch`: Symbol für Datenpunkte in einer Grafik (Kreis, Quadrat, Dreieck etc.) - `lty`/`lwd`: Darstellung (durchgezogen, gestrichelt, etc.) und Breite von Linien in einer Grafik - uvm. (siehe `?par`) --- ## Grafiken ### Mehrere Grafiken Für mehrere Grafiken in einem plot muss der Parameter `mfrow` genutzt werden: -- ```r par(mfrow = c(1, 2)) # Zwei Plots in einer Zeile (bzw. in zwei Spalten) plot(1:10, xlab = 'x', ylab = 'y') plot(1:10, xlab = 'x', ylab = 'y', type = 'l') ``` <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-14-1.png" width="600" style="display: block; margin: auto;" /> -- <br> Man beachte, dass dies eine *globale* Einstellung ist. Will man also wieder zum ursprünglichen Set-Up zurück, muss man mit `dev.off()` zurücksetzen. --- ## Grafiken ### Weitere Grafikenformen `hist()`, `boxplot()`, `barplot()` und `pie()`. -- #### Beispiel: ```r hist(mtcars$mpg, breaks = 5, freq = F, main = '', xlab = 'Miles per gallon') ``` <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-15-1.png" width="600" style="display: block; margin: auto;" /> --- ## Grafiken ### Funktionen, die grafische Elemente hinzufügen Es gibt auch Funktionen, die eine bestehende Grafik voraussetzen, bspw. -- - `lines()` - `points()` - `abline()` - `legend()` - `text()` - `arrows()` -- Bis auf `abline()` erwarten alle Funktionen X- und Y-Koordinaten. --- ## Grafiken ### Funktionen, die grafische Elemente hinzufügen &mdash; Beispiel ```r hist(rnorm(10000), freq = F, ylim = c(0, 0.5), xlab = '', main = 'Empirical vs. theoretical') x <- seq(-5, 5, 0.01) y <- dnorm(x) lines(x, y, col = 'red') ``` <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-16-1.png" width="650" style="display: block; margin: auto;" /> --- class: exercise_slide ## Grafiken ### Übungsaufgaben ```r set.seed(385) results <- rnorm(1000, mean = 100, sd = 15) ``` .font90[ 1. Kopieren Sie obigen Code und nehmen Sie an, dass dieser eine IQ-Testreihe mit 1000 Probanden simuliert. Zeichnen Sie ein Histogramm der Ergebnisse. Geben Sie Ihrem Plot anschließend eine passende Überschrift sowie passende Achsenbeschriftungen. Spezifizieren Sie darüber hinaus den Bereich von `x`- und `y`-Achse auf `\([40,160]\)` und `\([0,0.03]\)`. 2. Hinterlegen Sie dem Plot die dem IQ zugrundeliegende, theoretische Dichtefunktion, d.h. eine Normalverteilung mit `\(\mu=100\)` und `\(\sigma=15\)`. Wählen Sie als Zeichenfarbe Rot und machen Sie die einzuzeichnende Linie etwas breiter. 3. Zeichnen Sie einen Punkt in Form eines Dreiecks an das Maximum der theoretischen Dichte. Wählen Sie als Farbe Blau. 4. Kennzeichnen Sie sowohl das `\(0.025\)`- als auch das `\(0.975\)`-Quantil der theoretischen Dichte, in dem Sie Vertikalen an diesen Punkten einzeichnen. Wählen Sie als Farbe Grün. ] --- ## Lineare Regression ### Modell <!--Hauptaufgabe der Ökonometrie ist es, ?konomische Zusammenh?nge zu quantifizieren. Das Herzstück dabei bildet die Regressionsanalyse.--> .font80[ Einfaches (univariates) lineares Regressionsmodell: `$$Y_t=\beta_0+\beta_1 X_t+u_t,\quad t=1,\ldots,T,$$` wobei + `\(Y_t\)` - Regressand (zu erklärende Variable), + `\(X_t\)` - Regressor (erklärende Variable), + `\(u_t\)` - Fehlerterm. Mit KQ-Methode können wir die Koeffizientenschätzer `\(\widehat{\beta_0}\)` und `\(\widehat{\beta_1}\)` erhalten. ] -- .font80[ #### Beispiel: Im Datensatz **mtcars** sind reichlich Informationen zu verschiedenen Automodellen hinterlegt. Wir vermuten, dass die Reichweite eines Autos (*mpg*) davon abhängt, wie schwer es ist (*wt*). D.h. wir hätten folgendes Modell: `$$mpg_t=\beta_0 + \beta_1 wt_t + u_t$$` (Welche Richtung hat dieser Zusammenhangs vermutlich?) ] --- ## Lineare Regression ### Parameterschätzung </br> Einfache lineare Regression über `lm()` (= linear model) ```r model <- lm(mpg ~ wt, data = mtcars) ``` -- Die Ergebnisse werden in einer Liste gespeichert ```r names(model) ``` ``` ## [1] "coefficients" "residuals" "effects" "rank" ## [5] "fitted.values" "assign" "qr" "df.residual" ## [9] "xlevels" "call" "terms" "model" ``` --- ## Lineare Regression ### Regressionsoutput Einen kompakten Überblick über die wichtigsten Informationen erhält man mit `summary()` -- .font75[ ```r summary(model) ``` ``` ## ## Call: ## lm(formula = mpg ~ wt, data = mtcars) ## ## Residuals: ## Min 1Q Median 3Q Max ## -4.5432 -2.3647 -0.1252 1.4096 6.8727 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 37.2851 1.8776 19.858 < 2e-16 *** ## wt -5.3445 0.5591 -9.559 1.29e-10 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 3.046 on 30 degrees of freedom ## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446 ## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10 ``` ] --- ## Lineare Regression ### Punktewolke und Regressionsgerade ```r plot(mtcars$wt, mtcars$mpg, xlab = 'Weight', ylab = 'Miles per gallon') abline(model, col = 'red') ``` <img src="data:image/png;base64,#vorkurs_3_files/figure-html/unnamed-chunk-21-1.png" width="600" style="display: block; margin: auto;" /> --- ## Lineare Regression ### Vorhersage </br> Mit dem geschätzten Modell können wir Vorhersagen machen. Angenommen ein neues Auto mit einem Gewicht von `\(X_{new}=3\)` (in 1000 lbs) kommt auf den Markt. Was für eine Reichweite schätzen wir für dieses Auto? -- ```r new_weight <- data.frame(wt = 3) predict(model, newdata = new_weight) ``` ``` ## 1 ## 21.25171 ``` ??? was erwartet predict --- ## Lineare Regression ### Dummy-Variablen/Regression ohne Achsenabschnitt .font80[ Kategorische Variablen (hier Getriebeart *am* mit Automatik = `\(0\)`, manuell = `\(1\)`) können über die Funktion `factor()` als Regressor mit aufgenommen werden. ] -- .font80[ ```r model_wd <- lm(mpg ~ wt + factor(am), data = mtcars) model_wd$coefficients ``` ``` ## (Intercept) wt factor(am)1 ## 37.32155131 -5.35281145 -0.02361522 ``` ] -- .font80[ Will man eine Regression ohne Achsenabschnitt ( `\(\beta_0\)` ) durchführen, bedarf es einer `-1` (oder alternativ <code>+0</code> ) am Ende der Regressionsformel. ] -- .font80[ ```r model_wd_no_int <- lm(mpg ~ wt + factor(am) - 1, data = mtcars) model_wd_no_int$coefficients ``` ``` ## wt factor(am)0 factor(am)1 ## -5.352811 37.321551 37.297936 ``` ] --- class: exercise_slide ## Lineare Regression ### Übungsaufgaben .font90[ 1. Betrachten Sie im Folgenden den Datensatz `faithful`, der Daten zum [Old Faithful Geysir](https://de.wikipedia.org/wiki/Old_Faithful) im Yellowstone Nationalpark enthält. Sowohl die Dauer einer Eruption in Min. (`eruptions`) als auch die Wartezeit bis zur nächsten Eruption in Min. (`waiting`) sind als Variablen im Datensatz verfügbar. Unterstellen Sie nachfolgendes Regressionsmodell und schätzen Sie die entsprechenden Parameter und schätze die Parameter `\(\beta_0\)`, `\(\beta_1\)`. Zeichnen Sie anschließend eine geeignete Grafik und interpretieren Sie diese. `$$\text{waiting}_t=\beta_0+\beta_1 \text{eruptions}_t + u_t$$` 2. Verschaffen Sie sich mit `summary()` einen Überblick über ihr in 1 erhaltenes Ergebnis. Interpretieren Sie die geschätzten Koeffizienten und speichern Sie anschließend das `\(R^2\)` (*Multiple R-squared*) in `R2` ab. *Hinweise:* Schauen Sie sich die, beim Ausführen von `summary()`, ausgegebene Datenstruktur genauer an. 3. Angenommen Sie beobachten einen zusätzlichen Datenpunkt für die Dauer einer Eruption von `\(X_{new}=4\)`. Sagen Sie die entsprechende Wartezeit bis zur nächsten Eruption vorher. <!-- 4. Betrachten Sie nun den aus den Folien bekannten Datensatz **mtcars**. Regressieren Sie die Reichweite (*mpg*) auf die Getriebeart (*am* mit Automatik = 0, Manuell = 1) und lassen Sie den Achsenabschnitt weg. Interpretieren Sie ihr Ergebnis. --> ]